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Abstract 

In a recent work 10 , POULIN and one of us presented a quantum algorithm for preparing 
thermal Gibbs states of interacting quantum systems. This algorithm is based on Grover's tech- 
nique for quantum state engineering, and its running time is dominated by the factor D/ Zp, 
where D and Zp denote the dimension of the quantum system and its partition function at 
inverse temperature /?, respectively. 

We present here a modified algorithm and a more detailed analysis of the errors that arise 
due to imperfect simulation of Hamiltonian time evolutions and limited performance of phase 
estimation (finite accuracy and nonzero probability of failure). This modification together with 
the tighter analysis allows us to prove a better running time by the effect of these sources of 
error on the overall complexity. We think that the ideas underlying of our new analysis could 
also be used to prove a better performance of quantum Metropolis sampling by Temme et al. 

1 Introduction 

The ability to efficiently prepare thermal Gibbs states of arbitrary quantum systems at arbitrary 
temperatures on a quantum computer would lead to a multitude of applications in condensed 
matter, quantum chemistry and high energy physics |10 1 112 1 ^13], For example, we could estimate 
partition and correlation functions of fermionic and frustrated sytems. For these systems, the 
approach of first applying the "quantum-to classical map" [H] and then using the classical Monte 
Carlo method fails because the mapping does not conserve the positivity of statistical weights. 
We consider an arbitrary Hamiltonian H with spectral decomposition 

D 

H = Y,Eamm. (1) 

a=l 

The thermal Gibbs state of the system at inverse temperature f3 is given by 

Pp--=Y.-^\^a){M (2) 
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where := e denotes the partition function. 

The formal definition of preparing thermal Gibbs states is as follows: 

Problem (Thermalizing quantum states). Let H be a Hamiltonian, /3 an inverse temperature and 
€ G (0, 1) a parameter describing the desired accuracy. We consider the problem to prepare a state 
PfS that is e-close to the thermal Gibbs state p/s with respect to trace distanc^, i.e., 

\\p/3 - ppWtv < e. (3) 

We refer to the process of preparing such state as thermalizing the quantum system. We seek to 
determine efficient quantum circuits that realize such thermalizing process. 

We assume that the energies satisfy Ea £ [0, j]. If we initially only know that the spectrum 
of the Hamiltonian H is contained in the interval [£, u] , then the shifted and rescaled Hamiltonian 
4{H — iI)/{TT{u — l)) satisfies the condition of this assumption. The thermal Gibbs state is invariant 
under shifting of the spectrum. Thus, we have to rescale the inverse temperature by multiplying it 
by ?i — ^ when working with the new Hamiltonian. 

There are two types of quantum algorithms for preparing thermal Gibbs states. The first is 
a generalization of the Metropolis algorithm. The Metropolis algorithm can be applied to the 
special case of classical systems, i.e., systems whose Hamiltonian H = Yla ^a\o){a\ is diagonal in 
the computational basis. It offers great flexibility for constructing Markov chains whose limiting 
distributions are equal to the desired thermal Gibbs distributions. The number of times we have to 
apply the Markov chain scales like 1/6, where 6 is its spectral gap. Bounding the spectral gap from 
below for arbitrary systems and neighborhood structures is very difficult. However, it is possible 
to prove that the gap is sufficiently large for many practically relevant cases. 

Recently, Temme et al. presented an extention of the Metropolis algorithm to quantum sytems 
[12j . Their quantum Metropolis sampling makes it possible to implement quantum maps such 
that their fixpoints are approximately equal to the desired thermal Gibbs states. Analogously to 
classical case, the number of time we have to apply the quantum map depends on its spectral 
gap. The difficulty of bounding the gap from below remains for general systems and neighborhood 
structures. However, numerical experiments in [12] show that the gap scales like for the 
spin-chain Hamiltonian H = X^X^+i + Ifclfc+i + gZ^ on N spins. 

The second type of algorithm is due to Poulin and one of us [lOj. This algorithm behaves like 
a Las Vegas algorithm, i.e., it always produces a correct output p satisfying the requirements of the 
problem definition. The time it takes this algorithm to terminate is a random variable. However, 
we can bound the expected value. It is dominated by the factor \JD jZ^. This square root term 
occurs because this algorithm is based on an extension of Grover's state engineering technique. 

Once the Grover sampling has terminated we know that we have prepared a state that is close 
to the desired thermal Gibbs state. In contrast, we can only guarantee that quantum Metropolis 
sampling yields a good approximation if we have a lower bound on the spectral gap. But, of course, 
quantum Metropolis sampling has the potential to outperform the Grover sampling for certain 
quantum systems. 

We modify this Grover sampling and analyze the errors that arise due to imperfect simulation 
of Hamiltonian time evolutions and limited performance of phase estimation (finite accuracy and 
nonzero probability of failure) in more detail. This modification together with the tighter analysis 

^Recall that the trace distance is defined to be \tx\fX^X^ , where X = — pp. 
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allows us to prove a better running time. We show that the expressing the effect of these sources 
of error on the overall complexity is smaller than in the original algorithm. We also think that the 
ideas underlying of our new analysis could also be used to prove a better performance of the above 
quantum Metropolis sampling. 

This paper is organized as follows. In section [2] we present the structure of the algorithm. 
We identify three sources of errors that arise due to (i) imperfect simulation of Hamiltonian time 
evolution, (ii) limited precision of phase estimation, and (iii) non-zero failure probability of phase 
estimation. In section 13.11 and 13.21 we analyze how the complexity increases when we seek to keep 
the errors small. Finally, we make our conclusion in Section [H 



2 Quantum algorithm — idealized setting 

To better explain the intuition behind the quantum algorithm, we first ignore all sources of error. 
We assume that the unitary U = exp{2TTiH) can be implemented perfectly and efficiently. The 
eigenvalues Ea of H correspond to the eigenphases Ea of U, using the convention that the phase 
of e^'^*'^" is Ea- We assume that phase estimation (PE) makes it possible to perfectly resolve the 
eigenphases, i.e., there is an efficient quantum circuit mapping \ijja) <X> |00...0) onto \ijja) "X" \Ea) 
(this is the case as long as the energy Ea can be written as binary fractions). The realistic case is 
analyzed in detail in the following section. 

The algorithm prepares a purified Gibbs state of the form 

°-=^ ^ A B energy anc 

The states \(pa) form an orthonormal basis on the D-dimensional subsystem B. The \Ea) are 
computational basis states of the energy register, which consists of multiple qubits. These basis 
states encode the eigenvalues Ea of the eigenvectors iV'a) of H. The ancilla register consists of a 
single qubit. We obtain the thermal Gibbs state pp from by tracing out the subsystems B, 
energy, and anc (see eqn. (jl])) 

P;3=tr^(|/3)(/3|), (5) 

where we use A to denote the collection of the above three subsystems (the complement of A). 
The algorithms consists of the following steps: 
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Algorithm 1 Thermal Gibbs State Preparation at Inverse Temperature (3 

Input: Prepare the maximally entangled state \u) = W)^) o^i the subsystem AB. 

Step I: Run phase estimation of U on the A-part of \u). Write the eigenphase into the energy 

register. 

Step II: Apply the controlled rotation R = Y1e \ ^){^\ 'i^ Re where 

Re = 

The control is the energy register and the target is the ancilla qubit that is initialized in |0). 
Denote the resulting state by |^'). 

Step III: Use a variant of Grover to project onto the subspace in which the ancilla qubit is in 
|0). Denote the projector onto this subspace by Hq. The Grover iteration is given by 

G = (2|^)(^'| -/)(/- 2no). 

Output: The density matrix of final state |/3) by tracing out A. 





Let V be an arbitrary unitary. The maximally entangled state is invariant under the action 
oiViSiV, i.e., (y V)\h') = {u). Consequently, we can rewrite {u) as 



by setting V = J2a IV'a)(a| and \ipa) = V\a). In step I, we obtain the state 



In step II, we obtain the state 

I*) = ^ E (l^-) ® IV'-)) ® 1^-) ® [V^^\0) + Vl-e-^E^\l)) . (8) 



(9) 



' a 

Note that the desired purified Gibbs state \f3) is equal to 

nol^) 
l|no|^>|| ' 



where ||no|'I')|| = \J We apply the variant of Grover algorithm [2], which makes it possible to 
prepare \ j3) with an expected number of Grover iterations 0(l/||no|V') ||)- It is important that we 
do not need to know the overlap ||no|^)||. This shows that we obtain 



D 



in step III. 
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3 Quantum algorithm 



3.1 Analysis of simulation error 

The first source of error is the inabihty to implement U = exp(27rii7) perfectly for general H. 
Using techniques [TJ ITS] for simulating Hamiltonian time evolutions, we can only implement a 
unitary C/gim with \\U — Usha\\ < e^im- The resources grow inversely with the desired accuracy esim- 
To bound the error arising from imperfect simulation, we use the following result, which follows 
the discussion in [10', Appendix A]. 

Lemma 1. Let H he a Hamiltonian whose eigenvalues are contained in the interval [0,^]. Let 
U = exp{2TriH) and Usim be a unitary with \\U — UsimW ^ Csim- Then, there exists an effective 
Hamiltonian -ffsim such that C/sim = exp(27riffsim) o,nd \\H — -fTsimll < ^esim where k is a constant. 

Assume phase estimation could perfectly resolve the eigenphases of f/sim- Then, if we ran the 
algorithm using Us\m instead of U , then we would prepare the thermal state with respect to the 
effective Hamiltonian i^sim instead of H. Thus, it remains to determine how close the corresponding 
thermal states are close to each other with respect to trace norm. 

Lemma 2. Let H and //sim be as above. Then, the corresponding thermal states 

exp{-/3H) exp(-/3i7sim) 
and psim := —} . nrj — TV (11) 



tr(exp(-/3F)) tr(exp(-/3/7sim)) 
satisfy 

\\P - ftimlltr < 2 (12) 

provided that esim < e^/(8K/3). 

Proof. The fidelity of p and psim is given by 

F{p, Psim) = try'^A/^p^toV^ • (13) 
Using [3', Proposition 4] we bound the trace distance between p and psira as follows 



ll/O - /Osimlltr < Y 1 - -^(p,Psim)^ ■ (14) 

The analysis in \10\ Appendix C] shows that 

F{p,psim)>e-^^'^'^-, (15) 

and thus 

Hp - /Osimlltr <Vl- e-2/3«^sin. < ^2/3Kesim < |. (16) 

The rightmost inequality follows from 1 + x < e^ for all x G M. □ 

From now on, we measure the complexity in terms of how many times we have to invoke a 

controlled version of C/gim- If we wish to determine the complexity in terms of elementary gates, we 
have to look at the simulation technique more closely. 
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3.2 Analysis of Errors in Phase Estimation 



We now show how to prepare a state p such that \\psim — p||tr < e/2, implying that ||p — /o||tr < e as 
desired. We analyze the three phases of the algorithm. 



Phase I 

We need to run a special variant of phase estimation [8] of C/gim on \u). We briefly explain how it 
works. To avoid new definitions, we use iV'o) and Ea to refer to the eigenvectors and eigenphases 
of C/sirti) respectively. 

The usual phase estimation algorithm consists of the following steps [S]. The energy register 
consits of n = [log2(l/eprec)l qubits. We apply the Hadamard transform to each of the qubits of 
the energy register, the controlled- C/g^j^ gates (controlled by the j'th qubit of the energy register) 
on the ^-part of f, and the inverse quantum Fourier transform F'^ on the energy register. We 
measure the n qubits of the energy register in the computational basis and interpret the outcome 
b G [0,2*^ — 1] as the binary fraction E^, := 6/2", which is a very good estimate for Ea- More 
precisely, the probability of obtaining the estimate Ei, is given by 



1 sin2(7r2"(K - ^)) 
2^ sin2(vr(K-i5)) 



(17) 



We use \E}j) to denote the compuational basis state which encodes the energy value -Bfe. Let E^ 
denote the binary fractions that are closest to Ea, where we use the convention E" < Ea < Ea ■ 
It follows that the probability of obtaining Ea or Ea is greater or equal to ^ > |. Thus, the 
probability of failure, i.e., the probability of not obtaining one of the closest n-bit fractions, is less 
than |. 

To reduce the probability of failure to efaii, we repeat this quantum circuit k = [log2(l/efaii)] 
times, each time recording the estimate into a new energy register and adjoin a median register 
that consists of n qubits. This yields the state |T) 



1 



Eb 



i/ energy 



Eb 



where the amplitudes c 



Ea,Eiy 



satisfy |c 



Ea ,Ei, 



|2 = Pr(K,^6,)for 



J. /energy 



10. 



■0) median |0) anc ) 

(18) 



l,...,k. 



The median circuit determines the median of Eb-^ , ■ ■ ■ , Eb^. and writes it into the median register. 
Reordering the registers, we may write the resulting states as 



1^) = \^<^)\^<^) ® (Ca l^^>median «> \Pa 



' energy 



+ l'^a)median(g)energy®") |0)j 



(19) 



where 



the states |//^) are supported only on the states |£'f,J 
Eb^ , • • • , Eb^. is equal to E^ , and 

the states |^a) are orthogonal to \Ea) (8> |/^a )• 



) such that the median of 
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It follows from the analysis in [51 Hj that the amplitudes satisfy 

l-efaii< |c+P + |c-|2 < 1, 0< |||ea)f <efaii. (20) 

This means that the probability of the median not being one of the closest binary fractions E'^ to 
Ea is less than or equal to efaii- The advantage of combining phase estimation with the powering 
techinque for approximation algorithms is that we only need to invoke a controlled version of C/gim 

r(l/eprec)log(l/efail)l (21) 

instead of 0((l/eprcc)(l/efaii)) when using phase estimation alone [9]. 

To keep the notation simple, we use \E^) to denote the tensor product \E^) \fJ,^). Using this 
convention, we write the state after step I (phase estimation) as 

^ a 

Phase II 

The i?-operation is controlled by the energy value contained in the median register. After step II, 
|^>) evolves to the state 



1^) = 4^Y.\'f'a)\^'^)(^ct\Et)0{V^^\O) + ^/l^7^\l))+ (23) 

^ a 

^ V ' 

^^m\iPa)(^Rm<s)\o)) . (24) 

^ a 

' V ' 

10 

As a consequence of the property {(,\ip) = 0, we have 

^ Y^iCal^a) < efaii and 1 - efaii < m)f < 1. (25) 



2 ^ 

D 



Phase III 

Let 1/3) be the state obtained by applying Grover's algorithm to |^'), i.e. 



|no|*)|| 



(26) 



Let p be the reduced density operator of over A. We need to bound ||no|^')|| from below to 
obtain an upper bound on the expected number of Grover iterations. We also need to show that p 
is close to Psim- This is done in the following lemma. 
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Lemma 3. Let p = tr^(|/3)(/3|). This density operator has the form 



^ tr^(no|V;)(^|no) ^ tr^(no|0(e|no) 
(§|no|§) (^'inol^') 



that satisfies 



II ~|| .|| tr^(no(|V;)(V;|)no) „ , „ tr^(no(|0(g|)no) „ .e e_e 

Psim — P tr S Psim , ~ , ~ , tr + , = , ^ , tr S 7 + 7 — - , i^oj 

(^'IIIol^') (^'IIIol^') 4 4 2 

provided that efaii = e~^ and eprec = e/(32/3). 

Proof. Observe that the off-diagonal terms no(|V')('?| + IO(^l)no in (no|^)(^|no) vanish when we 
trace over A. 

Set := tr(no|^)(^'|no) and define the operator 

a := tr^(no|V)(V|no) = ^ ^(Ic^ pe-'^^"^ + \c-\\-P''^ )\^a)m • (29) 

a 

We can express = tr(c7) + tr(no|0(e|no). Since tr(no|0(e|no) = (eiHolO < M? < efaii, by ([25]) 
we can bound and cr's trace as follows 

(1 - efaii) e-^^p-^ < tr(a) < e^^^"^ , (30) 

(1 - efaii) ^ e-'^^f'-^ < iV < ef,n + ^ e^^i'-^ . (31) 
Because G [0, -j], we can bound the ratio Z[f3)/D as follows 

1 > Z(/3)/Z) > e^^. (32) 
By choosing the lower bound on and the upper bound on ||tr^(no(|0(^|)no)||tr5 we obtain 

tr^(no(|0(g|)no) ^^f^^:^^^l^^^^^^2Re,...^^_ 

^ (l-efail)t"(l-^fail)- 



tr^(no(|V.)(V.|)no) _ " e-P^^ tr(a) e 



The last inequality is obtained because e is small and e^ <1 + 2x for x G [0, 1]. The term 

d 

is still satisfied even when examining the following two extreme cases 

(I) Lower bound on A^ and upper bound on tr (cr): : 3 1 , (34) 

^ ^ (1 - efaii)e-/^^p- ^ ^ 

n _ efaii) e~^^p''°'=) 

(II) Upper bound on A^ and lower bound on tr((7): 1 . (35) 

^ efaii + e^^p-<^ 
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We know that 



Because 1 + 2x > for Vx € (0, 1), we derive 

e^^p-<^ 1 + 1 e 

7 ^ — -R 1 < — - 1< T- (37) 

(1 - efaii)e-/^^p-- - 1 - efaii "4 

In the second case because 1 + x < for Vx G M and 1 < D/Z{f3) < e^, we have 

^ _ (1 - efaiQe-^^"-) ^ ^ _ e-/^^p-- ^ ^ _ 1 - ]| ^ e 

+ e^^p--- ~ £2 + e/^^p-'^ ~ £2 + 1-4 ^ ^ 

for small e. □ 

4 Conclusion 

Theorem. Let H be a Hamiltonian, (3 an inverse temperature and e G (0, 1) a parameter describing 
the desired accuracy. Let f/sim &e the Hamiltonian simulation such that \\U — UsiraW ^ Csim where 
U = exp{2TTiH). Our algorithm prepares a state pp that is e-close to the thermal Gibbs state pjj, 
i.e., 

11/5/3 - ^/3||tr < e, (39) 

provided that egim < e^/(8/3K), eprec = e/(32/3) and efau = e'^e"^. The complexity of our algorithm 
scales like 

ofJ^-(log- + /3)) (40) 



e e 

in terms of the number of invocations of the controlled-U sim operation. 

Proof. The requirements for esimj eprec and efaii are immediate by Lemma [Hand Lemma By (pT]) 
the cost for performing one Grover iteration scales as 

O (log i + /?)). (41) 

The number of Grover iterations ^ is determined by 0{j^^^) = 0{^J ^) when using the lower 
bound of tr((T) in pO]) . 

5 Acknowledgments 

P. W. and C. C. gratefully acknowledge the support of NSF grants CCF-0726771 and CCF-0746600. 



9 



References 

[1] D. Berry, G. Ahokas, R. Cleve and B. Sanders, Efficient Quantum Algorithms for Simulating 
Sparse Hamiltonians, Communications in Mathematical Physics, vol. 270, pp. 359-371, 2007. 

[2] M. Boyer, G. Brassard, P. Hoyer and Alain Tapp, Tight Bounds on Quantum Searching, 
Fortschritte Der Physik, vol. 46(4-5), pp. 493 - 505, 1998. 

[3] C. Fuchs and J. Graaf, Cryptographic Distinguishability Measures for Quantum Mechanical 
States, IEEE Transactions on Information Theory, vol. 45, issue 4, pp. 1216-1227, 1999. 

[4] M. Jerrum, L. Valiant and V. Vazirani, Random Generation of Combinatorial Structures from 
a Uniform Distribution, Theoretical Computer Science, vol. 43, issue 2-3, pp. 169-188, 1986. 

[5] P. Kaye, R. Lafiamme, and M. Mosca, An Introduction to Quantum Computing, Cambridge 
University Press, 2007. 

[6] S. Lloyd, Universal Quantum Simulators, Science, vol. 273. no. 5278, pp. 1073 - 1078, 1996. 

[7] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of State 
Calculations by Fast Computing Machines, J. Chem. Phys., vol. 21, pp. 1087-1092, 1953. 

[8] D. Nagaj, P. Wocjan, Y. Zhang, Fast QMA Amplification, QIC vol. 9 no. 11&12 pp. 1053-1068, 
2009. 

[9] M. Nielsen and I. Chuang, Quantum Computation and Quantum Information, Cambridge 
University Press, 2000. 

[10] D. Poulin and P. Wocian, Thermalizing Quantum Systems and Evaluating Partition Functions 
with a Quantum Computer, Physical Review Letter, vol. 103, pp. 220502, 2009. 

[11] M. Suzuki, cd. Quantum Monte Carlo Methods in Equilibrium and Nonequilibrium Systems, 
vol. 74 of Springer Series in Solid-State Science, Springer, 1988. 

[12] K. Temme, T. Osborne, K. VoUbrccht, D. Poulin and F. VerstraeteK, Quantum Metropolis 
Sampling, arXiv: abs/0911.3635, 2009. 

[13] B. Terhal and D. DiVincenzo, Problem of equilibration and the Computation of Correlation 
Functions on a Quantum Computer, Physical Review A, vol. 61, pp. 022301, 2000. 

[14] P. Wocjan, C. Chiang, D. Nagaj and A. Abeyesinghe, A Quantum Algorithm for Approximating 
Partition Functions, Physical Review A, vol. 80, pp. 022340, 2009. 

[15] C. Zalka, Proc. R. Soc. London, Ser. A, Simulating Quantum Systems on a Quantum Computer, 
vol. 454, no. 1969, pp. 313 - 322, 1998. 



10 



